Time-dependent density-matrix functional theory for biexcitonic phenomena 
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We formulate a time-dependent density-matrix functional theory (TDDMFT) approach for higher- 
order correlation effects like biexcitons in optical processes in solids based on the reduced two- 
particle density-matrix formalism within the normal orbital representation. A TDDMFT version of 
the Schrodinger equation for biexcitons in terms of one- and two-body reduced density matrices is 
derived, which leads to finite biexcitonic binding energies already with an adiabatic approximation. 
Biexcitonic binding energies for several bulk semiconductors are calculated using a contact biexciton 
model. 
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I. INTRODUCTION 



The theoretical description of ultrafast processes in 
modern electronic devices is an important problem of 
contemporary condensed matter physics*^ The entan- 
gled role of the fluctuation and correlation effects, espe- 
cially in low dimensions, makes such an examination a 
challenge. In particular, it is not easy to reproduce cor- 
rect excitonic and biexcitonic features in the optical ab- 
sorption spectra of materials. Besides fundamental inter- 
est such as in 4- wave mixing jiaa excitonic and biexcitonic 
effects have a variety of practical applications, such as 
optoelectronic devices' entangled photon sources' and 
quantum computing^ 

The standard approaches, based on the semicon- 
ductor Bloch equations (SBEs) 8 and nonequilibrium 
Green's function techniques'"^ cannot be easily ap- 
plied to study higher-order correlation effects in strongly 
nonequilibrium situation, because this requires many- 
particle correlation functions that depend on many time 
arguments. 11 Approaching these problems with time- 
dependent density-functional theory (TDDFT)^ looks 
promising due to its formal simplicity and the fact that 
it in principle includes correlation effects exactly; how- 
ever, this usually requires going beyond the standard 
LDA-GGA approximations^ Excitonic effects have been 
studied with TDDFT in several ways, including the time- 
dependent optimized effective potential approach^ and 
the combination with the Bethe-Salpeter method^ Un- 
fortunately, these approaches also become very tedious 
in the strongly nonequilibrium case. 

Recently, we proposed an alternative approach for ul- 
trafast excitonic effects based on the single-particle den- 
sity matrix and a TDDFT version of the SBEsJ^^ 
We showed that the effective electron-hole attraction is 
defined by matrix elements of the exchange-correlation 
(XC) kernel / xc with respect to the valence and con- 
duction band Kohn-Sham single-particle wave functions. 
Experimentally observed lowest exciton binding energies 
can be reproduced in a simple way using the local and 



long-range XC kernels 

f^(r,r') = -A S(r-r') (1) 

and 

f^') = -^Ty ( 2 ) 

where Aq and a can be viewed as adjustable parame- 
ters. The ALDA / xc leads to a too weak electron-hole 
attraction to produce bound excitons. 

In the case of biexcitons, which are correlated dou- 
ble electronic excitations, the problem is much more 
complicated. One reason is that multiple excitations 
in TDDFT require nonadiabatic XC functionals, and so 
far there are no simple approximations available. An- 
other reason is that at first sight it is not clear how 
to represent biexcitonic wave functions in Kohn-Sham 
TDDFT. In this paper, we formulate and test an alterna- 
tive TDDMFT approach for biexcitons based on the nat- 
ural orbital (NO) representation for the stationary elec- 
tron eigenf unctions where the multiparticle excited 
states are naturally related to the higher-order density 
matrix elements. 

II. TDDMFT FORMALISM 

The standard TDDFT single-particle Hamiltonian is 

h(r,t) = ~ +V(r,t) + V B [n](r,t) + V«[n](r,t) ■ (3) 

Here, V(r, t) — U nuc i(r) + Kxt( r , t) is the static potential 
of the nuclei plus the time-dependent external perturbing 
potential. We consider a homogeneous external electric 
field in dipole approximation, T4 x t(r, t) = —rE(t), which 
implies that the characteristic field frequency is much 
larger than the level spacing^ Vn(r,i) is the Hartree 
potential, and V xc (r,t) is the time-dependent xc poten- 
tial, which are both functionals of the time-dependent 
single-particle density n(r,t). The Hartree potential is 
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not very important for the description of excitions, but 
the xc potential is crucial, since it accounts for the effec- 
tive electron-hole interaction.— 

In general, Vx C [n](r,i) has a memory, i.e., it depends 
on densities at previous times t' < t. The resulting xc 
kernel / xc , defined as 



/ xc (r,r',w) = 



[ > 5n(r>,t') 



(4) 



»o(r) 



therefore has in general a frequency dependence. An ex- 
plicitly frequency-dependent / xc is required for describing 
double excitations with linear-response TDDFT* 21 How- 
ever, to date there are only few approximations for / xc 
available, and none of them is particularly suited for the 
biexcitonic properties in solids we have in mind. For this 
reason, we choose a slightly different approach. 

To describe the properties of doubly-excited iV-elec- 
tron systems, one can consider the one- and two-electron 
density matrices, defined asi£~— 



j(xi , x' x , t) = N J dx 2 J dxz... J dxN 

x ty(x 1 ,X2,... ,x N ,t)^f*(x' 1 ,x 2 ,... ,x N ,t), (5) 

T(xi, X2, x[, x' 2 , t) = N(N — 1) / dxj, / dx^... / dx 



x ^(xi,x 2 ,-,x N ,t)^*(x' 1 ,x' 2 ,...,x N ,t), (6) 

where 'J is the many-body wave function and Xi = (rj, s,) 
denotes the space coordinate and spin index. 7(0:1, x[, t) 
describes the single-particle properties of the system, 
such as the charge density n(x, t) = ^(x, x, t). Moreover, 
all ground-state quantities, including Tq{x\,X2,x' 1i x' 2 ), 
can in principle be obtained from 70(21, 2^), since there 
is one-to-one correspondence between 70 {x\ , x[ ) and the 
ground state many-body wave function ^0 (the DMFT 
generalization of the Hohenberg-Kohn theorem 2 ^) . 

Let us now restrict the discussion to two-electron sys- 
tems and derive equations of motion for the density ma- 
trices. We consider the following effective two-electron 
Hamiltonian: 

H(r u v 2 ,t) = h* d {v u t) + h ad (r 2) t) +Mn 2 ]( ri ,r 2 ,t),(7) 

where ft ad is the TDDFT Hamiltonian © using an 
adiabatic approximation for V xc [n](r,t), which leads to 
a frequency-independent xc kernel /^(r, r'). In this 
way, excitons can still be described, since the frequency- 
dependence of / xc is not essential for the electron-hole 
interaction. However, biexcitons (which are correlated 
two-particle excitations) cannot be captured in the adi- 
abatic approximation. To make up for this, we intro- 
duce, in a somewhat ad-hoc manner, an effective two- 
particle interaction u;[?i 2 ](ri, r 2 , t) which we define as 
a functional of the two-particle density n 2 (ri,r 2 ,i) — 
v P*(ri, r 2 , t)^{vi, r 2 , t). In this way, dynamical screen- 
ing effects can in principle be accounted for, as done in 
standard many-body perturbation theory. 



In the following we express the two-electron wave- 
function in terms of the NOs Xfc( r )- I n the singlet case 
one obtains *(r,r',i) = J2k,i Cki(t)xk(r)xi(r'), where 
Cki(t) is a symmetric matrix, and k, I are the appropriate 
quantum numbers (band index, momentum, spin etc). 
We shall use this matrix in general non-diagonal form for 
physical insight on the nature of the excitations. Since 
the density matrices r y(xi,x[,t) and T(xi,x 2 ,x' 1 ,x' 2 ,t) 
are defined by the two-electron wave function, they can 
be expressed in terms of the matrix elements Cki(t): 



7(ari,a4,t) = ^7fci(*)XfcO c Ox*( a:; i)> ( 8 ) 

T(x 1 ,X 2 ,X / 1 ,X 2 ,t) = ^ ^klmn(t)Xk(xi)xi(x2) 
klmn 

*Xm(zi)Xn04), (9) 



where jki(t) = ^Y, m C k m {t)C^{t) and T k i mn (t) = 
2Cfci(t)C^„(i). We will soon see that in the two-band 
approximation the excitonic wave function is propor- 
tional to 7k" k 2 (t) , and the biexcitonic one to r^^k k4 (t) , 
where c and v stand for the conduction and valence 
bands, and is the corresponding electron and hole mo- 
mentum. From now on, we shall use superscripts to de- 
note band indices. 

The equation of motion for 7fe;(i) and Tkimn{t) can 
be obtained via the equation for Cki(t). From the time- 
dependent two-electron Schrodinger equation one finds: 



% dCkdt)_ = j2(h kr (t)c rl (t) + C kr (t)h rl (t)) 



dt 



^ WklrsC rs (t) 



(10) 



with the initial condition Cki(t = 0) = SkiCk and the 
matrix elements 



h kr (t) = I dr X Ur)h ad (r,t)xr(r) 



(11) 

Wki mn {t) = J dr i J d r 2Xfc(riM( r 2)Mn 2 ](ri,r 2 ,i) 
XXm(ri)x„(r 2 ) . (12) 

Here and in the following, it is implied that each spatial 
integration is divided by the unit cell volume. Equa- 
tion (|10p is nonlinear, since the matrix elements, which 
depend on the electron density, are functions of C k i(t). 
From the definitions of jki (t) and Tkimn it) one then ob- 
tains the following equations for the one- and two-particle 
matrix elements: 



£^2M = y^^hkrjrl - Jkrhrl) 



dt 



(^krsm w msrl ^krsmW m srl) , (13) 
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Tklrs)- (14) 



An important feature of Eqns. (fT3|) and (fl4[) is the fact 
that they are closed, i.e. one does not need to truncate an 
infinite hierarchy of equations for higher-order density- 
matrix elements. However, keep in mind that this prop- 
erty is only valid for two-level (two-band) systems. In 
the single electron (w = 0) linearized diagonal approxi- 
mation for two bands, one obtains the TDDFT-Wannier 
equation for the exciton eigenenergies and eigenfunctions 
from Eq.(7)J£ 

KA« = E [K'+q - 4<) <W + F kk ,] 7&, iq , (15) 



where k is the electron momentum, q is the sum of the 
electron and hole momenta (the exciton momentum) and 

F kk <-2 J dvj dr' X : k (r) Xwk (r')/ X c(r,r')x; k '(r')Xck'(r') 

(16) 

are the matrix elements for the effective electron-hole at- 
traction. 



III. TWO-LEVEL MODEL FOR BIEXCITONS 

The possibility to obtain a biexcitonic state with the 
TDDMFT formalism can be already shown for a two-level 
model with the energy levels E\ and E 2 > Ex and the 
HOMO-LUMO gap E g = E 2 - E x . From Eqs. (JTSJ) and 
(|14p in the lowest (second) order approximation (by keep- 
ing only those matrix elements which contain no more 
than two indices "2"), one obtains the following system 
of equations after carrying out a Fourier transformation 
into the frequency domain: 



(u-E g - fV 1 - G 1 r 2211 = 

(u) - 2E g - G 2 )r 2211 



_ F 2 7 21 = Q) 
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(17) 
F 2 = 



where F 1 = W2121 — Willi + W2112 + ti- 2 \ + w : 
2h 2 i + W2212 + W2221 + W21H Gl = ^12 + W2122 - W1112 
<4m, G 2 



W2222 - wini 



w 22 n, and 



hl b = / dld2 X l{l) Xb {l) 



5Kc(l) 



w 



abed 



dld2d3 X :(l)x* b (2) 



8n{2) 
5w{\,2 



X2(2)xU2), (18) 



6n(S) 



x Xc(l)x«i(2)xa(3)x*(3), 
w r abcd = [ dld2d3d4x* a (l)xt(2) 



(19) 



6w(l,2) 

x Xc(l)Xd(2)X2(3)x 2 (4)xt(3)x^4). 




V — exo 



FIG. 1: (Color online) Excitonic bound states in a two-level 
system with single-particle energies E\ and E-2- E cxc and 
Fbiexc are the binding energies of excitons and biexcitons, 
respectively. A biexcitonic state can be thought of as aris- 
ing from a two-step process: First, two excitons are created, 
which then combine to form a biexciton whose energy is less 
than the energy of the two individual excitons. 



All matrix elements are evaluated at the initial (non- 
perturbed) densities, and we use the shorthand notation 
1 for ri, 2 for r 2 , etc. The solutions of the system of 
equations (ITT1) can be easily found: 



3E a + F 1 



G 2 1 
~ ± 2 



{E g ~F 1 + G 2 ) 2 +AG 1 F 2 . 

(21) 

From the general solution, one can discuss several lim- 
iting cases. In particular, for no correlations (F 1 = F 2 — 
G 1 = G 2 = 0) one gets a trivial solution with one and 
two free excited electrons: u>i = E g , u) 2 = 2E g . In 
the absence of a density-dependent two-electron poten- 
tial (w[n 2 ](ri, ri) = G 2 = 0), one finds an excitonic state 
with energy uj\ = E g — F 1 = E g — E cxc < E g . However, 
the second root oj 2 cannot be lower than 2E g — 2E 0XC , 
i.e. the linear adiabatic approximation does not produce 
a biexciton, as expected. On the other hand, for nonzero 
w[n 2 ](ri,ri), when G 2 < and IF 1 ] < \G 2 
obtain a biexcitonic level with binding energy G 

Therefore, in order to obtain a biexcitonic state in pure 
TDDFT [with no w[n 2 ](ri, ri)] in the adiabatic approx- 
imation, one needs to go to the nonlinear regime and 
consider the next order terms (~ 7r) in the equations. 
To show the possibility of a biexcitonic solution one as- 
sumes that the initial state includes a long-living exci- 
ton, 7r — > 7T, where 7 is the averaged "excitonic func- 
tion" (see also Ref. where a possibility to obtain 
double excitations in the adiabatic TDDFT was consid- 
ered). Then, the second (biexcitonic) equation (ITTl) ac- 
quires an additional term (at G 2 =0), which results in 
the eigenenergy 2E g + 2(h 22 — hj^j that corresponds to 

- ^11)7 <0 an d 
Possible excita- 



one can 
2 2jF u 



the biexcitonic solution in the case (h 2 
2E g + 2\(hl 2 -hJ 1 )j\>2\F 1 \=2E c 



(20) tions in the two-level case are illustrated in Fig. [TJ 
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IV. THE TWO-BAND CASE 



The generalization to the many-electron case is straightforward. The equation for the biexciton energies (with zero 
center-of-mass momentum) in the two-band approximation has the following form: 



= 



■ d c v v ' 

*^ — £ k+q — £ k' + £ k + £ k'+q 



1 k+q,k',k,k'+q — 2-^ U k+q,l<';k+q,k',k,k'+q i k+q,k\k,k'+q ' ^ > 



k' 



where 



C ----- — n cv I A 1 - --- -l- A 2 - --- 14- 4 3 - --- 

"^k+q.k'ik+q.k'.k^k'+q — ^k+q.k' ^k';k+q,k',k,k'+q T ^k+q^k'lk+q^'.k.k'+qj ^ ^k+q.k'ik+q.k'.k.k'+q 

+ (k + qok'), (23) 
^k; kl ,k 2 ,k3,k 4 = ydlrf2d3x: k (l)x: k (l)ffi(l,2,3)xck 1 (2)Xck 2 (3)x,k 3 (2)*x: k4 (3), (24) 

^.k^kaaa, = / dlrf2d3d4^ k (l)Xck'(2)Xck(l)^k'(2)52(l,2,3 ) 4)x ckl (3)x ck2 (4)x !) k 3 (3)*x;k 4 (4), (25) 



^k,k';k 1> k 2 ,k 3 ,k 4 = 2 A k,k';k 1 ,k 2 ,k3,k 4 [Xck(ri) -+Xvu(r2% (26) 

and gi(r,ri,r2) = sn^rfll) anc ^ 9z{ T U r 2, r 3, r -i) = tnllllll) are two-particle density kernels. Similar to the excitonic 
case, Eq. (|22[) is the momentum representation version of the Schrodinger equation for two electrons and two holes,— 
where the matrix elements G k+q k /. k+q ^ k k / +q correspond to an integral inter-particle (in general, four-body) inter- 
action. 

To solve Eq. (|!22"|) . we expand the biexcitonic function in terms of the complete set of the excitonic functions 7™ k 
with eigenenergies £n,q (n is the number of the bound state), which can be found from the solution of Eq. (fTB"j) . and 
antisymmetrize it with respect to interchange of holes and electrons, in order to satisfy the Pauli principle. Then the 
biexcitonic functions can be expressed in the following form: 

-pccVu'i \ ^ v ^.v' l± -j- v v' t± (07\ 

k+q,k',k,k'+q / j m,k-fq,q 'm.k'+q,— q y nm,q ' 7n,k',k' — k 7m,k+q,k— k fU nm,k f - k ' v / 

where ± correspond to two possible states of biexcitons, singlet (— ) and triplet (+), with respect to two-electron 
spins.— Thus the problem is reduced to finding the matrix elements that enter into Eq. (J22J and using the orthogo- 
nality of the excitonic eigenfunctions, one finds the equation for the biexcitonic eigenvectors and the corresponding 
eigenenergies, similar to Eq. (|15p : 



where 



^ [i u ~ E nq - E m q) <W <5mm"W _ H nm,n> m' ,qq> b n'm',d' - °> ( 28 ) 



n^m^q' 



^Tirrt,fi J m J ,qq' (["^ ~F ^] )nm,n'm\ qq' ' ( 2 *^) 

Snm,n'm',qq' — 7n,k+q,q7m,k+q',-q7n',k+q',q'7m',k+q.-q' ! (30) 

k 

and 

W ± = -Sjnm'Sgg, ^ 7n,k+q.q- F k,k'7n'.k+q'.q' ~ <W<W ^ 7m,k,-q^k,k'7m'.k',-q 
k,k' k.k' 



-)- \ ^ E 1 * J ^f* n'* v v' , v' '* v v' 

/ j k.k 7 I 7n,k / +q,q7m,k+q', — q7n',k+q',q / 7m / ,k+q / , — q ' 7m,k',- q7n,k— q'+q,q7n / ,k,q / 7m',k— q'+q,— q' 

k,k' 

^ 7n,k+q.q7m,k',-q [^k+q,k';k+q',k',k,k'+q' T Gk+q.k' M> M+q> MM> +q>] 7^', k +q',q'7m', k ',-q'- (^1) 



k,k',k,k' 

r 



Equation (|28[) formally resembles equation (|15[) for ex- citons, with the exciton eigenenergies used instead of 
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TABLE I: Excitonic and biexcitonic binding energies E^ xc and 
^biexc £ Qr gome semiconductors (in meV), calculated with the 
TDDMFT formalism. The parameters Aq, Ai and A2 which 
determine the model interaction kernels {TJ, (|32p and (|33p 
have been adjusted so that the calculations reproduce the 
experimental exciton binding energies. 





A 


77'CXC 


C A 1 /Q 2 


A 2 /& 


zribicxc 


ZnO 


290 


60 


1.82 


101 


15 


CdS 


308 


28 


0.022 


0.64 


5.7 


CuCl 


20.7 


190 


1.97 


37.2 


32 


CuBr 


20.9 


110 


2.0 


11.5 


25 



the bare single-electron energies. The effective exciton- 
exciton attraction H^ m n , m , , is defined by the matrix 
elements of the kernels gx(ri, r 2 , r 3 ) and <jf 2 (ri,r 2 , r 3 , r 4 ). 

To test the formalism, we consider the simple case of a 
single biexciton in one-exciton level approximation, n — 
m = 1 (see also Fig. [TJ . We have obtained the solution of 
this equation using simple model kernels: the local kernel 
(P) to generate the excitonic states, and the following 
local two-particle kernels for biexcitons: 

gi ocal (r, n, r 2 ) = -C A 1 8{v - v 1 )6{v - r 2 ) (32) 

[which includes the averaged element Co of the one- 
electron excited density- matrix component C cv , see 
Eqns. p5 ]) -([2l) j) ] and 

g^\v,v',v 1 ,v 2 ) = -A 2 5{v-v')5{v-v 1 )8{v-v 2 ). (33) 

The kernels (I3"2l and (j3"3")l can be viewed as constituting 
a "contact biexciton" model, in analogy with the contact 
exciton model defined by the xc kernel ([T]).— 

Results for the electron eigenenergies and eigenfunc- 
tions of several semiconductors were obtained by using 
the VASP 4.6 code2£ with GGA-PAW potentials and a 
350 eV energy cutoff. We approximated the NO func- 
tions by the corresponding Kohn-Sham single-particle 
wave functions, which can be considered as a good ap- 
proximation when the correlations are not too strong. 
We find that with the effective local kernels / x ° cal , g l ° cai 
and .g 2 ocal one can reproduce the experimental biexcitonic 
binding energies with a suitable choice of the parameters 
Aq, Ai and A 2 (see Table I). 

Let us briefly discuss how one can in principle find 
the nonadiabatic kernel / xc (w) from the adiabatic poten- 
tial iu[n](ri,r 2 ) which produces biexcitonic states. For 
this, one can use an approach similar to the one pro- 
posed by Maitra et al. for double excitations in the two- 
electron case<21 Namely, the expression for / xc (a;) can be 
obtained by expanding the excited two-electron excited 
wave function \&(ri, r 2 , i) in terms of quasi-degenerate 
wave functions of the single-particle excited and biex- 
citonic states, and by comparing the eigenenergy equa- 
tion [for the Hamiltonian (J7J, which includes w[n](ri, r 2 )] 



with the corresponding TDDFT Casida equation. The 
single-electron excited states have to include all states 
which are close to the biexcitonic one and are well sep- 
arated from the other states. A detailed formulation for 
such a general case will be reported elsewhere. 



V. CONCLUSION 

In this paper we have formulated a TDDMFT ap- 
proach to study biexcitonic effects. We have derived 
the TDDMFT version of the Schrodinger equation for 
biexcitons in terms of the two-particle density-matrix el- 
ements in the two-band approximation. We have solved 
this equation in the case of several semiconductors by 
using phenomenological two-electron interaction kernels, 
thereby defining a contact biexciton model. With this 
model one can reproduce the lowest biexcitonic bind- 
ing energies by using proper kernel parameters. To ob- 
tain biexcitonic states within the single-particle TDDFT 
approach, one would either need to use a frequency- 
dependent XC kernel or consider the nonlinear regime. 
Generalization for the case of bound states with larger 
number of particles is in principle straightforward. 

There are several advantages of this simplified formal- 
ism for biexcitons comparing to other approaches: 1) 
it can be adapted for use in the real-time domain in a 
straightforward manner; 2) physical transparency of the 
method, in particular the effective TDDFT electron-hole 
and exciton-exciton interactions are directly related to 
the interaction kernels, which may allow one to make sim- 
ple estimations of the possibility to produce bound states 
with given TDDFT kernels; 3) in many cases it may al- 
low one to construct the non-adiabatic Kohn-Sham XC 
kernel from the phenomenological adiabatic two-particle 
density kernel, which may shed some light on the general 
requirements on the Kohn-Sham f xc (ui) necessary to pro- 
duce biexcitons and other higher-order coupled states. 

Examination of ultrafast processes and higher-order 
correlation effects, including the excitonic and biexcitonic 
transport, in semiconductor nanostructures and organic 
molecules is underway. 
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